† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant Nos. 10732010 and 10972010).
The three-dimensional premixed H2-O2 detonation propagation in rectangular ducts is simulated using an in-house parallel detonation code based on the second-order space–time conservation element and solution element (CE/SE) scheme. The simulation reproduces three typical cellular structures by setting appropriate cross-sectional size and initial perturbation in square tubes. As the cross-sectional size decreases, critical cellular structures transforming the rectangular or diagonal mode into the spinning mode are obtained and discussed in the perspective of phase variation as well as decreasing of triple point lines. Furthermore, multiple cellular structures are observed through examples with typical aspect ratios. Utilizing the visualization of detailed three-dimensional structures, their formation mechanism is further analyzed.
Detonation is a supersonic flow with chemical reactions. Since discovered in an experiment of flame propagation, it has been studied for over one hundred years. Nowadays the research on mechanism of detonation ignition and propagation, especially on structural analysis of cellular patterns, is still popular due to its importance in construction blasting, production security, military facilities,[1] and so on.
Limited by experimental equipment and technique, it is difficult to capture the details of detonation structures. As a result, smoked foil technique[2] and high-speed schlieren technique[3] were developed in 1960s to observe the detonation structure indirectly. White et al.[4] was the first to make a detailed experimental research on three-dimensional gas detonation. Strehlow[5–7] gave a well-documented description of detonation structure. In his papers, the structure of detonation cells had been described as an interaction of Mach configurations. Hanana et al.[8] experimentally obtained clear patterns recorded on the wall of square ducts. Through analysis, they considered that there were two kinds of stable cellular structures in three-dimensional cases: the diagonal mode and the rectangular mode. Lin et al.[9] improved the experimental set to study the critical structure in the transition process of the deflagration to detonation (DDT) under different boundary conditions.
Compared with experiments, numerical methods have advantages in providing global distribution of any variable at any time. A key problem of simulating the detonation is how to describe the chemical reaction in numerical method. Taki et al.[10] used a two-step kinetics model to analyze the structure of transverse wave. Oran et al.[11] proposed a detailed reaction model examining the cellular pattern of two-dimensional oxyhydrogen detonation. Besides, one-step kinetics model[12] and other simplified models[13] were also put forward in previous work. Sichel et al.[14] developed Taki’s model by introducing a parameter which reflected changes in gas composition, and promoted the computational precision a lot. In the early study, two-dimensional (2D) simulations were usually used to study the detonation cellular structures and there were thousands of literatures available. However, detonation was an intrinsically three-dimensional process. With the rapid development of computational fluid dynamics (CFD), three-dimensional simulations gradually became feasible. Williams et al.[15] adopted the one-step reaction model and numerically modeled three-dimensional structures in self-sustaining detonations. Their simulation started from a 2D periodic solution, and added a one-dimensional perturbation perpendicular to the motion of the initial 2D transverse waves. Tsuboi et al.[16,17] used a Harten–Yee non-monotone upstream-centered scheme to simulate three-dimensional hydrogen/air detonation and obtained the rectangular mode and diagonal mode. Deledicque et al.[18] also presented a detailed numerical study of three-dimensional structures in gaseous detonations using the one-step reaction model and a parallelized, unsplit, shock-capturing algorithm. Dou et al.[19,20] and Wang et al.[21,22] used the one-step kinetics model and the fifth order weighted essentially non oscillatory (WENO) scheme to simulate detonation propagation in square tubes. They found that the detonation wave in small size of square tubes could propagate spirally just as in round tubes, which was named as spinning mode. Additionally, Weng,[23] Ivanov,[24] Cai,[25] and Huang[26,27] also did a lot of simulation work with three-dimensional rectangular ducts and obtained similar results. However, to the best of the authors’ knowledge, most of previous work focusing on three-dimensional problems required a relatively high-order scheme and cost a plenty of computing resources.
In this paper, we employed the in-house parallel detonation code based on an improved three-dimensional space–time conservation element and solution element (CE/SE) scheme with second-order accuracy[28] to further study the three-dimensional cellular structures of detonation propagating in rectangular ducts. For the first time, we attained the cellular patterns with critical cross-sectional size as well as different aspect ratios.
Detonation is an extreme chemical reactive flow involving strong shocks. Therefore, to simulate the detonation process accurately, the following two aspects are crucial: (i) an appropriate chemical reaction model and (ii) a good shock capturing scheme. Currently, detailed chemical reaction models[11] are most close to the real situation, since elementary reactions actually occur during the detonation. However, it requires too many computational resources even for 2D cases, not to mention three-dimensional computation. To achieve an acceptable compromise between the solution accuracy and the computational time, we choose an improved two-step kinetic model,[14] in which the complex chemical reactions are simplified to an induction reaction and an exothermic reaction. Neglecting the influences of thermal conduction, viscosity, and mass diffusion, the detonation propagation can be described in three-dimensional Cartesian coordinate system by the reactive Euler equation
For shock capturing, we employ the CE/SE scheme,[29] which is a unique computational framework for solving hyperbolic equations. This numerical method takes time and space both as equal mathematical dimensions. Conservation elements (CEs) and solution elements (SEs) are defined appropriately to make local and global conservation law be strictly guaranteed. CE is a space–time control volume on which the conservation law is implemented. Meanwhile, SE is used to solve the fluxes involved in the discrete conservation law. The dissipation of the scheme can be efficiently controlled, because it is constructed based on a non-dissipative scheme.[29] As a result, it not only has a good shock capturing capacity, but also can be used to simulate acoustic waves.[30] In order to solve problems with complex boundary conditions as well as with a lower expense of computational resource, unstructured[31] and essentially conservative adaptive[32] definition of CEs can be also employed. But these improved schemes will lose some degree of global conservativeness. Moreover, dimension splitting method is not required for multi-dimensional cases, resulting in a genuinely multi-dimensional scheme.[28,33] These properties make CE/SE method be very efficient for detonation simulations. Up to now, it has already been successfully used for simulating gaseous detonation[34,35] and liquid-fueled two-phase detonation.[36,37]
Notably, the characteristic time scales for flow and chemical reactions may have several orders’ difference. To save CPU time, a decoupling method is used to handle the stiff source terms. That is to say, we first solve the flow field without considering chemical reactions and then integrate the source terms using sub-time step. According to the Gauss divergence theorem, the governing equation (
![](cpb_25_11_114702/cpb160693ieqn1.gif)
![](cpb_25_11_114702/cpb160693ieqn2.gif)
![](cpb_25_11_114702/cpb160693ieqn3.gif)
![](cpb_25_11_114702/cpb160693ieqn4.gif)
![](cpb_25_11_114702/cpb160693ieqn5.gif)
![](cpb_25_11_114702/cpb160693ieqn6.gif)
![](cpb_25_11_114702/cpb160693ieqn7.gif)
Numerous factors, including the size and aspect ratio of cross-section, the initial perturbation, the reaction model and so on, exert a coupling effect on detonation structures in a rectangular duct. We focus on the effect of the size and aspect ratio of cross-section (AR) herein. Although the cellular structures in square tubes have been comprehensively studied before, we carry out similar simulations both as a code validation and as baselines for comparison. The settings of 14 total cases are listed in Table 1. There are five things to note. (i) At least 50 grid points are guaranteed for one cellular width to ensure the accuracy of simulation results. We have tried beforehand the simulation with grid resolution of 10, 25, and 50 points per millimeter (pts./mm). After comparison of the cellular patterns, we concluded that 25 pts./mm is enough for examples with cross-sectional size greater than 2 mm, while it had better be 50 pts./mm when cross-sectional size equals to 1 mm. (ii) It needs to lengthen the computational domain to 160 mm when the cross-sectional size reduced to 2 mm or 1 mm, because detonation structure at first 80 mm in these cases is not always stable. (iii) All the initial perturbations are only given at first 30 time steps to trigger transverse waves. (iv) Each example takes 8 CPU cores in parallel running. (v) In all cases, high pressure is initially implemented on the left 5% part of the computational domain (P = 45 atm and T = 298 K) and the rest part is filled with H2–O2 mixtures of normal temperature and pressure (P = 1 atm and T = 298 K). Additionally, reflective boundary condition is performed for all the side walls, the left wall of CPU 1 and the right wall of CPU 8.
![]() | Table 1. Simulation parameters in different cases. . |
In the cases with square cross-sections, the diagonal mode and the rectangular mode are observed in Cases 1 and 2 with relatively large cross-section size, and the spinning mode is obtained in Cases 5 and 6 with very small size, as shown in Figs.
![]() | Fig. 1. Pressure isosurfaces and motion of TLs during about one period in (a) Case 1, (b) Case 2, (b) Case 3, (d) the middle domain of Case 4, (e) Case 5, and (f) Case 6. |
Logically, critical structures appear in Cases 3 and 4 with medium cross-section size. In Case 3 with initially diagonal perturbation, the elementary cells shown in Fig.
Furthermore, we extract the data of MPHs along the central and angular lines in x direction and make the diagram displayed in Fig.
Generally speaking, when the size of cross-section is big enough, it will form diagonal or rectangular pattern, depending on types of initial perturbation. Otherwise, when the size of cross-section is very small, the detonation structures will finally convert into spinning patterns, however we change the other conditions. Additionally, critical structures also exist while the diagonal or rectangular mode is transformed into the spinning mode. It is the process with accumulation of phase difference and decreasing of TLs. Comparing with the rectangular mode, the diagonal one is more difficult to be transformed into the spinning mode, since TLs in this condition originally move along the diagonal line, different from spinning structures.
In different propagation modes, the AR makes different influences on cellular structures. Firstly, we assess the effect with the cell size equal to 4 mm. Cases 1, 7, and 9 with AR ranging from 1 to 3 are carried out to study multiple cellular structures in the diagonal mode. When the length in one transverse direction of square tubes is stretched integer times, detonation cells generally duplicate the same number, as in Case 9. But if carefully examining the cellular structure of Case 7 in Fig.
![]() | Fig. 4. Cellular patterns on the lateral walls with typical ARs in (a) Cases 1, 7, and 9, as well as in (b) Cases 2, 8, and 10. |
![]() | Fig. 5. Distribution of MPHs with motion of TLs on y–z slices in Case 7 with (a) simulation results, and (b) one or (c) another kind of imaginary structures, as well as (d) in Case 8 and (e) Case 10. |
Similarly, Cases 2, 8, and 10 are used to study multiple cellular structures in the rectangular mode. For all the cases, the number of cells on sides of the ducts does not change like Case 7, as shown in Fig.
While simulating with the cell size of 1 mm, the effect of ARs differs a little from examples above. A case study of initially diagonal perturbation and AR = 2 (Case 11) is shown in Fig.
![]() | Fig. 6. Cellular patterns on the surface in Case 11, including (a) y–z slices and (b) expansion of lateral walls. |
On the other hand, results in Case 12 with initially rectangular perturbation have something different. Observing the cellular patterns in Fig.
![]() | Fig. 7. Cellular pattern on x–y and x–z sides in Case 12 with x ranging from 0.075 m to 0.089 m (a), and from 0.127 m to 0.141 m (b). |
![]() | Fig. 8. Two kinds of ∞ shaped trajectory of the cellular patterns on y–z slices in Case 12. (a) One direction of ∞ shaped trajectory. (b) The other direction of ∞ shaped trajectory. |
The most remarkable difference between structures in Figs.
As expected, if the AR becomes big enough, the detonation structure will degenerate to the 2D situation. Take simulation results in Case 14 as an example, in which we set AR = 4:1 with initially rectangular perturbation, as shown in Fig.
An in-house parallel detonation code based on the CE/SE scheme with second-order Taylor expansion as well as two-step kinetic model is used to simulate the three-dimensional detonation in rectangular ducts with premixed H2–O2 atmosphere. We chiefly obtain three well-understood types of propagation mode as well as critical structures through typical examples in square tubes, and then examine the influence of AR on detonation cellular structures in the rectangular, diagonal and spinning mode. Ultimately, the simulation reveals that: (i) the process of the rectangular or diagonal mode transforming into the spinning mode occurs with phase variation as well as decreasing of TLs; (ii) multiple cellular structures will generally consist of elementary cells in square tubes with certain phase difference or periodic inequality, or degenerate to the 2D situation when assigned a high value of AR. Based on our results, we conclude that quantitative change in size or aspect ratio of cross-section will finally make qualitative differences on the three-dimensional detonation cellular structures.
1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 |